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Using large-scale dynamical cluster quantum Monte Carlo simulations, we explore the unconven- 
tional superconductivity in the hole-doped Hubbard model on the triangular lattice. Due to the 
interplay of electronic correlations, geometric frustration, and Fermi surface topology, we find a 
doubly degenerate singlet pairing state at an interaction strength close to the bare bandwidth. Such 
an unconventional superconducting state is mediated by antiferromagnetic spin fluctuations along 
the Y-K direction, where the Fermi surface is nested. An exact decomposition of the irreducible 
particle-particle vertex further confirms the dominant component of the effective pairing interaction 
comes from the spin channel. Our findings provide support for chiral d + id superconductivity in 
water- intercalated sodium cobaltates Na^CoCb • 2/H2O, as well as insight into the superconducting 
phases of the organic compounds ft-(ET) 2 X and Pd(dmit)2. 
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Introduction.- Since the discovery of the Cu-based 
high temperature superconductors, the search for new 
unconventional superconductors is among the central 
topics in condensed-matter physics [1]. The water- 
intercalated sodium cobaltates Na^CoC^-^/H^O [2|-y] and 
two families of organic charge-transfer salts ft-(ET)2X 
and Pd(dmit)2 [5l-tlQl[ are of particular interest. The un- 
derlying structure of these layered materials is the ge- 
ometrically frustrated triangular lattice. The competi- 
tion between electronic correlations and geometric frus- 
tration yields novel phenomena. For example, the most 
frustrated members of the tt-(ET)2X and Pd(dmit)2 fam- 
ilies are believed to host quantum spin liquid states |ll| . 
and the recently discovered 5K superconducting phase 
in Na^CoC^.^/H^O might be a chiral state which breaks 
parity and time reversal giving rise to interesting edge 
modes that can carry quantized particle and spin cur- 
rents dMii. 

The layered triangular lattice compound Na^CoC^ • 
2/H 2 has a superconducting dome for x ~ 0.3, y ~ 1.3 at 
T C ^5K [2|-y|. Due to intercalation, its electronic struc- 
ture is effectively two-dimensional. A very rich phase 
diagram has been mapped out for a range of Na concen- 
trations [4]; however, the nature of the superconducting 
phase has remained poorly understood. Recent measure- 
ments on high quality single crystals [15( show that the 
spin contributions to the Knight shift decreases below T c 
along the a and c axes, supporting the notion that the 
Cooper pairs are formed in a spin-singlet state. The tem- 
perature and doping dependence of the Knight shift and 
the relaxation rate above T r provide evidence of antifer- 



romagnetic correlations [15L llq 



There are a number of theoretical proposals for the 
unconventional superconductivity in the cobaltates. The 
underlying triangular lattice allows a doubly degenerate 
E2 representation of the superconducting order param- 
eter with d x 2_ y 2 and d xy degenerate states [l|, LL3|, Il7j . 



raising the exciting possibility of a time-reversal symme- 
try breaking chiral d x i_ y i ±id xy superconductor. Earlier 
studies of the cobaltates draw analogy to the cuprates 
and employed either a phenomenological RVB theory 
[18], or a slave boson mean- field approach [19| to pro- 
vide signatures of a spin-singlet d-\-id pairing state. Re- 
cent studies of the sodium cobaltates using a Gutzwiller 
projection supplemented by symmetry arguments [12], as 
well as the multi-orbital functional renormalization group 
[14j . reveal a rich phase diagram with an anisotropic d-\-id 
phase and a possible topological quantum phase transi- 
tion through a nodal superconducting state. However, 
prior approaches suffer either from their mean-field na- 
ture, or their incapability of capturing correlation effects 
in the strong coupling regime. Hence, there is an urgent 
need of unbiased studies, where the interplay of strong 
electronic correlations and geometric frustration can be 
treated in a non-perturbed fashion. 

The simplest model that captures the essential physics 
of the cobaltates is the single-band Hubbard model on 
a triangular lattice. In this Letter, we explore the low- 
energy properties of this model by large-scale dynamical 
cluster quantum Monte Carlo simulations [20]. We fo- 
cus on the different superconducting instabilities in the 
hole-doped side of the phase diagram. To the best of 
our knowledge, this is the first study of the hole-doped 
Hubbard model on the triangular lattice exploring the 
pairing symmetries on different cluster sizes. Clusters 
up to size N c = 12 allow a greater momentum resolu- 
tion and higher quality data on the spectral function, 
self energy, and different superconducting susceptibili- 
ties. Therefore, we obtain an unambiguous signature 
of an unconventional doubly-degenerate superconducting 
state in the strong to intermediate coupling region. By 
explicitly comparing the pairing susceptibility in the s-, 



d X y-waxe singlet channels and the /-wave triplet 
! and d xy components are 
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channel, we find that the cL 



most divergent and extrapolate to the same T c within 
our numerical accuracy. We identify that the pairing is 
mediated by strong spin fluctuations along the antifer- 
romagnetically (AF) ordered wavevector on the r to K 
direction. The Fermi surface (FS) is nested along this AF 
wavevector, but the system only orders at half- filling in 
the Heisenberg limit. An exact decomposition of the ir- 
reducible particle-particle vertex furthermore reveals the 
dominant part in the effective pairing interaction comes 
from the spin channel. 

Formalism.- The Hamiltonian of the system is H = 
Eka e k c L c ka + U y £ i n it n ii , where cj^c^) is the cre- 
ation (annihilation) operator for electrons with momen- 
tum k and spin a, n^ = c^Cia- is the number operator, 
and the bare dispersion is given by e£ = — 2tcos(k x ) — 
4tcos(y/3k y /2)cos(k x /2) with t being the hopping am- 
plitude between nearest neighbor sites, and U the on-site 
Coulomb repulsion. 

We investigate one and two-particle properties of 
the model using the dynamical cluster approximation 
(DCA) [2l| with weak-coupling continuous time quan- 
tum Monte Carlo (CTQMC) [22] as the cluster solver. 
The DCA maps the original lattice onto a periodic clus- 
ter of size N c embedded in a self-consistently determined 
host. Spatial correlations inside a cluster are treated ex- 
plicitly while those at longer length scales are described 
at the mean-field level. In this work we choose clusters 
of sizes N c = 4, 6, 8 and 12. We study inverse tempera- 
tures up to fit = 16.5. We obtain the cluster self-energy 
E(K,cj) via the maximum entropy method [23| (MEM) 
applied directly to the Matsubara- frequency self energies 
calculated by the DCA-CTQMC 0, [Hj. We then in- 
terpolate the E(K,cj) to obtain the lattice self energy, 
E(k, lj), and lattice spectral function, A(k,uj). 

To obtain various susceptibilities, x(T), we extract the 
irreducible vertex function T via the Bethe-Salpeter equa- 
tion from the two-particle Green function measured on 
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n = 1.0 
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the cluster, then employing x(^) = 

1 - Txo 
is the bare susceptibility constructed from the dressed 

one-particle lattice Green function. The superconduct- 
ing pairing susceptibilities are obtained from the particle- 
particle channel, and the charge and spin susceptibilities 
are obtained from the particle-hole channel. We further 
separate the pairing susceptibilities explicitly into spin 
singlet and triplet channels, where in the singlet channel 
we project the x(T) pairing onto s-, d x 2_ y 2-, and d^-wave, 
and in the triplet channel we project it onto the /-wave 
channel, with the corresponding form factors |lL Il4j. 

To explore the pairing mechanism we decompose the 
particle-particle pairing vertex T into the fully irreducible 
vertex A, the charge (S = 0) particle-hole contribution, 
<1> C , and the spin (5 = 1) particle- hole contribution, <£ s , 
through the parquet equation, r = A + <1> C + <l> s [26]. We 
furthermore project the previous expression using differ- 
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FIG. 1: (Color online) Cluster spin susceptibility for N c = 12, 
interaction strength U — 8.5£, temperature T — O.lt and 
different fillings, n = 0.667, 0.8, and 1. 
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where each term is the projected component of the corre- 
sponding term in the parquet equation [27]. In this way, 
we are able to distinguish which component contributes 
the most to the effective pairing interaction. 

Results. - Fig. [1] displays the cluster spin susceptibility 
at different fillings, n = 0.667, 0.8, and 1, and coupling 
U = 8.5t. The data are obtained from DCA-CTQMC 
simulations with cluster size N c = 12, and we interpolate 
the cluster susceptibility into the entire Brillouin zone 
(BZ). At very high hole-doping, n = 0.667, the suscepti- 
bility is mostly flat. As the filling increases, n = 0.8, the 
spin susceptibility develops six bumps at the K points. 
When n = 1, the bumps become more pronounced. The 
vector connecting r to K is the antiferromagnetically or- 
dered wave vector (Qaf) in the Heisenberg limit of the 
half-filled model. The cluster spin susceptibility demon- 
strates that the antiferromagnetic fluctuations become 
stronger as the filling moves towards n = 1. The pair- 
ing of electrons may be mediated by these fluctuations 

Fig.[2]shows the Fermi surface (FS) at the same fillings 
used in the previous figure. Fig. [2] (a) corresponds to the 
non-interacting limit. At n = 0.667, 0.8 and 1 the FS 
is close to a perfect circle. The van Hove singularity in 
the non-interacting band structure is present at n = 1.5 
with saddle points at M. One- loop RG calculations [30( 
show the FS in the hole-doped side is stable against weak 
Coulomb interactions. However, under strong interac- 
tion, the FS begins to deform. Fig. [2] (b) displays the 
FS at n = 0.667, U = 8.5t, which is slightly deformed 
towards a hexagon. The red arrow corresponds to Qaf, 
while the pink arrow is this vector shifting its center to 
T and rotating it by 60°. For n = 0.667 the pink ar- 
row is longer than the diameter of the FS so there is no 
nesting effect, and we do not observe superconductivity 
at this filling. In Fig. [2] (c), n = 0.8, the FS is more 
deformed, the Qaf now connects significant sections of 
the FS, and as illustrated in Fig. [TJ the AF fluctuations 




FIG. 2: (Color online) (a) First Brillouin zone, the symmet- 
ric path r — M — K — r and the non-interacting Fermi sur- 
face at different band fillings, (b), (c), (d) Spectral function 
A(k,u = 0) on the Fermi surface for N c = 12 DCA-CTQMC 
simulations with U = 8.5£, T = O.lt, and n = 0.667 in (b), 
n = 0.8 in (c), and n = 1 in (d). Red arrow is the AF ordered 
wave vector (Qaf) and the pink arrow is after shifting its 
center to F and rotating it by 60°. 



are stronger. The nesting effect and the strong AF fluc- 
tuations together give rise to diverging pairing suscepti- 
bilities at filling n = 0.8 and 0.9, as discussed below. At 
half-filling, n = 1, the FS is further deformed towards 
a hexagon, but the spectral weight become less coherent 
due to the strong interaction. Interestingly, the nesting 
vector now is shorter than the diameter of the FS. Hence, 
even through the AF fluctuations are the strongest here, 
electrons on the FS are hard to pair by Qaf, the sys- 
tem is rather subject to a Mott transition, whose novel 
features are beyond the scope of this paper. 

Fig. [3] displays the inverse pairing susceptibility as a 
function of temperature, l/Xpairing(T), at filling n = 0.9, 
U = 8.5t, and N c = 6. Here we explicitly project the lat- 
tice pairing susceptibility in the s-, d x 2_ y 2- 1 and d xy -wdive 
singlet channels and the /-wave triplet channel by using 
the appropriate form factors. Fig. [3] shows that the two 
singlet d-wave components are the most divergent ones. 
Within our numerical resolution their 1/Xpairing extrap- 
olate to zero at the same superconducting transition tem- 
perature, T c . This implies that the superconducting 
order parameter is doubly degenerate with components 
d x 2_ y 2 and d xy . Based on symmetry arguments any lin- 
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FIG. 3: (Color online) Inverse pairing susceptibility, 
l/Xpairing, for N c = 6, U — 8.5£ and n — 0.9. The singlet 
s-wave and triplet /-wave do not diverge, whereas the singlet 
pairing channels with d x 2_ y 2- and d xy -wave symmetry show 
a divergencency at the same T c . Note that we have multiplied 
by a factor of 7 the s-wave pairing susceptibility in order to 
use the same vertical scale. Inset, the superconducting transi- 
tion temperature T c as a function of doping. FL and SC label 
the Fermi liquid and superconducting regions, respectively. 



ear combination of both d-wave components is possible 
below T c . However, both Ginzburg-Landau and BCS- 
type mean-field approaches favor superconducting phases 
that break the time-reversal symmetry for singlet multi- 
component superconductors |3lL |32| , 



such as the d + id 



singlet pairing state predicted in graphene [131 . Il7j and 
the cobaltates 12l. Il4j. Therefore, our findings support 
the possibility of a chiral d + id superconducting phase in 
the hole-doped triangular Hubbard model. The inset of 
Fig. [3] shows the phase diagram for different doping con- 
centrations based on N c = 6 DCA-CTQMC simulations. 
T c becomes finite for doping larger than n = 0.7 due to 
the onset of FS nesting and strong AF correlations, and 
increases as n aproaches 1, reflecting that the AF fluctu- 
ations become stronger towards half-filling. However, the 
nature of the ground state at half-filling is still unclear 
due to a worsened minus-sign problem in our simulations. 
To shine light on the dominant contribution to the 
pairing interaction, we use the parquet equations to de- 
compose the irreducible particle-particle vertex function, 
and project each term onto its d x 2_ y 2 and d xy compo- 
nents. The results are presented in Fig. |4] for a DCA- 
CTQMC simulation with cluster size N c = 12, U = 8.5t 
and filling n = 0.9. The left, right panels correspond 
to the d xy: d x 2_ y 2 projection of the parquet equations, 
respectively. In both cases, the dominant contribution 
to the effective pairing interactions Vd and Vd 2 _ 2 i s 
from the magnetic, spin 5 = 1, particle- hole channel, 
Vf and Vf . In fact, we also find that the pair- 

d xy « x 2_ y 2 ' ^ 

ing interaction, V dx j d 2 _ 2 (k — k') is peaked at momen- 
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FIG. 4: (Color online) Left: d xy projected contributions to 
the pairing vertex Vd xy , from the fully irreducible vertex V dx , 

charge V dx and spin v£ x cross channels versus T at n = 0.9, 
C/ = 8.5£. Right: the d x 2_ y 2-wave projection of the same 
quantities. In both cases, the contribution to the pairing in- 
teraction from spin channel is clearly dominant. 



turn transfer |k — k'| = |Qaf|- The vertex decomposi- 
tion confirms that this peak comes from the spin channel 
V Z Id 2 2 (Qaf) (not shown). Note that both N c = 6 
and 12 size clusters have the cluster joints connected by 
Qaf From the BCS gap equation 



=4i> sc ( k - k ') 



A, 



2E(k') 



.^(^1), 



(2) 



where E(k) = \f^ + A£, we infer that if the super- 
conducting pairing interaction V sc '(k — k') is peaked 
at Qaf, the order parameters Ak which correspond to 
d x 2_ y 2 , d xy and /-waves are equally favored in the N c = 6 
and 12 clusters. Our results suggest that d xy and d x 2_ y 2 
singlet pairing are favored over the /-wave triplet pair- 
ing, probably because /-wave has a more complex nodal 
structure than the two d- waves [l|. 

Conclusion.- Using large-scale dynamical cluster 
quantum Monte Carlo simulations, we find a doubly de- 
generate singlet pairing state at interaction strength close 
to the bare bandwidth and filling larger than n = 0.7 in 
the hole-doped Hubbard model on the triangular lattice. 
Our findings suppport the presence of a chiral d + id 
singlet superconducting phase in this model. The pair- 
ing mechanism comes from antiferromagnetic spin fluc- 
tuations at the magnetic order wavevector nesting the 
deformed FS. A decomposition of the vertex further con- 
firms that the spin channel contributes the most to the 
effective pairing interaction. 
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